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We address how the nature of hnearly dispersing edge states of two dimensional (2D) topological 
insulators evolves with increasing electron-electron correlation engendered by a Hubbard like on-site 
repulsion U in finite ribbons of two models of topological band insulators. Using an inhomogeneous 
cluster slave rotor mean-field method developed here, we show that electronic correlations drive the 
topologically nontrivial phase into a Mott insulating phase via two different routes. In a synchronous 
transition, the entire ribbon attains a Mott insulating state at one critical U that depends weakly 
on the width of the ribbon. In the second, asynchronous route, Mott localization first occurs on the 
edge layers at a smaller critical value of electronic interaction which then propagates into the bulk 
as U is further increased until all layers of the ribbon become Mott localized. We show that the 
kind of Mott transition that takes place is determined by certain properties of the linearly dispersing 
edge states which characterize the topological resilience to Mott localization. 
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Topological insulators (TI) are a new quantum phase 
of matter distinguished by a nontrivial topology of their 
electronic state[l-7]. The TI phase has been predicted 
and discovered in numerous compounds starting from the 
two-dimensional quantum spin hall (QSH) systems[3, 4] 
to various three-dimensional materials[8~12]. While the 
search is on for new materials with exotic topological 
character[13, 14], the phenomenon presents opportunities 
to unravel new physics. 

While the essential features of topological insulators 
(TI) are manifestations of one-electron physics, there is 
intense interest to explore the effect of electronic corre- 
lation and disorder on the topological phase[15-21]. In 
TIs, strong spin-orbit coupling (SOC) leads to a time- 
reversal invariant band structure with a bulk charge gap 
and gapless edge modes which show Dirac like linear dis- 
persion with a characteristic velocity[7]. The gapless 
edge states are immune to weak perturbations (disor- 
der/interactions) that are time reversal symmetric. A 
natural question that arises is regarding the fate of 
the Dirac dispersion with increasing electronic interac- 
tions/correlation. The system is expected to evolve to 
a Mott insulating state with increasing local repulsion 
characterized by a scale U . What is the mechanisim of 
such a Mott transition, in particular, does the edge mode 
velocity renormalize to zero, in an analogous fashion as 
the Brinkman-Rice[22] mechanisim of a diverging effec- 
tive mass in the Hubbard model? Previous works[15-20] 
which addressed similar issues revealed that the gapless 
surface states do survive weak to moderate correlation 
while under stronger correlations the TI phase evolves 
into various insulating phases e.g., a spin liquid phase 
with gapped charge spectrum but gapless surface spinon 
excitation or a bulk antiferromagnetic insulator etc. 

Of particular interest is the role of the local repulsion U 
on the nature of electronic states in finite systems of TIs 
with terminating edges (we focus on 2D in this paper), 
such as, for example, a ribbon which is long along the 



x-direction but of finite width L along the y-direction. 
The question of how interaction affects differentially the 
gapless edge states and the gapped bulk states existing 
simultaneously in such system of TIs with finite bound- 
aries is unexplored hitherto. Such a study requires a 
detailed treatment that captures inhomogeneous nature 
of the electronic state due to the lack of translational 
symmetry perpendicular to the boundaries. This, along 
with the necessity to tackle strong interactions makes the 
problem difficult, even prohibitively expensive, to treat 
with accurate methods such as exact numerical diago- 
nalization [17]. Here we study interaction effects in TI 
ribbons by employing an inhomogeneous cluster slave ro- 
tor mean-field (SRMF) method[23-25] which is known 
to provide a correct qualitative description of the Mott 
physics [15, 25] while allowing for the treatment of large 
system sizes. The effect of local interaction is captured by 
introducing a Hubbard like on-site repulsion {U) into two 
well known models of 2D TIs - the Kane-Mele (KM)[1, 2] 
and Bernevig-Hughes-Zhang (BHZ)[3] models. Our focus 
is on the dynamics in the charge sector and the concomi- 
tant Mott localization engendered by increasing U , fol- 
lowing the spirit of ref . [22] . 

The highlights of the study are the following. With in- 
creasing U , a finite ribbon of TI attains the Mott insulat- 
ing state in two distinct ways. First, in the synchronous 
transition, charge fluctuations vanish simultaneously on 
all layers (with different y-coordinates) of the ribbon at a 
critical Uc that depends on the ribbon width L. The sec- 
ond route to the Mott state is via an asynchronous tran- 
sition. Here, the charge fluctuations at the edge layers 
vanish at a critical Ue which is independent of L. With 
further increase in U , successive layers become Mott lo- 
calized until at an L dependent Uc the entire ribbon is 
Mott insulating. A remarkable feature here is the si- 
multaneous realization of two phases in the ribbon - the 
Mott insulating (topologically trivial) state in outer lay- 
ers and topologically non-trivial state in inner layers for 
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Ue < U < Uc- This results in gapless modes in the bulk 
layers that are at the boundary between the Mott insu- 
lating and topologically non-trivial regions. We discuss 
the physics underlying these results later in the paper. 
TI models: We briefly describe the two models for topo- 
logical insulators that we study here. The KM model[l, 2] 
describes electrons hopping on a honeycomb lattice, 
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where t is the nearest neighbor, hopping ampli- 

tude and A is spin-orbit coupling strength. de- 
note next nearest neighbor sites, is 2:-Pauli spin ma- 
trix and = ±1 depending on the orientation of the 
hop[l]. For A 9^ 0, the resulting dispersion has non-trivial 
topology. The BHZ model [3] defined on a square lattice 
with four spin-orbit coupled atomic orbitals, e.g., |st)i 
\P t) = \{Pv + iPx) t), \s i), and \p \,) = \{py - ip^) |) is. 
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where a,/? = s,p and a =^,i. ta is the energy of the 
spin-orbit coupled orbital \aa) and 5 denote a nearest 
neighbor vector. The hopping matrix elements tsa,ap in 
the |scr), \pa) basis is given by (see also [26]), 
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where tss, tgp, tpp are overlap integrals and a takes values 
-1-1 (—1) for spin f (i). We set tss = tpp — t and tsp/ \/2 = 
v/2, and define eq such that eg = — (eq ~ 4t) and tp — 
(eo — 4i). Thus we are left with three parameters for 
the model - t^ v and cq. The model shows inverted band 
structure and hence the topological phase for < eq < 8i. 
Interactions: We introduce electron interaction via an 
on site Hubbard repulsion, 
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to the two models in Eq. (1) and (2) where Ui is the 
electron number operator. For the KM model = 



cl^Ci^i while for the BHZ it is rii 
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Inhomogeneous SRMF Formulation: We now out- 
line the inhomogeneous cluster SRMF formulation de- 
veloped here to study the above two interacting mod- 
els for lattices with ribbon geometry. In the usual slave 
rotor (SR) formulation [23-25], the electron operator is 
expressed as a product of a fermionic spinon (/) which 
carries the spin and a bosonic rotor (9) which carries the 



charge. One writes q^^. = fl^^e ' subject to the con 



straints, J2ao 



= 1 and J2a <a = Eaa 



terms of the auxiliary operators is then mean-field de- 
coupled as follows. Positing the ground state of Hsr to 
be the direct product of the spinon and rotor ground 
states, 1^) — l^*/) l^'e), one defines effective spinon 
and rotor Hamiltonians as = 7^^^ [vPe) and 
n" = (4-/1 H^^l^-/), respectively. For the BHZ model, 
this gives 

'^B=Y^i<^o,~ l^f)fLaf^o,a~ t^fLaf^+SfSa (5) 
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where chemical potentials are introduced to control the 
mean particle density. The effective hopping parame- 



ters tf and t^ are given by t{g^^^^ = tsa,apBi^i+s and 
where Bi 
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and x1a,jp = (*/l//aa/j73<T I*/) (analogous expressions 
are used for the KM model). The resulting spinon and 
rotor Hamiltonians are solved self consistently to ob- 
tain the ground state. A key point to be noted here 
is that the Bs and xs are bond dependent in a rib- 
bon, i. e., they are inhomogeneous. The spinon Hamil- 
tonian is quadratic albeit with inhomogeneous renor- 
malized hopping parameters and hence can be diago- 
nalized numerically[27] The rotor problem, however, is 
non-quadratic and a full numerical solution is still pro- 
hibitive, and additional approximations such as cluster 
mean field method have to be adopted. In conventional 
cluster mean-field method [25], one considers a small clus- 
ter of sites where the rotor Hamiltonian is treated exactly. 
Terms connecting the sites inside the cluster to sites out- 
side (bath) are decoupled by using, <I>i = (e*^') where 
$i is a site-dependent mean-field parameter also to be 
obtained self-consistently. This quantity plays an impor- 
tant role in that it signifies charge fluctuation at the site 
i. A nonzero $i implies charge fluctuation whereas van- 
ishing of implies a "local" Mott insulating phase. 
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FIG. 1. Inhomogeneous cluster mean-field method. Each 
unique bond in the strip which is treated as a distinct cluster 
is shown enclosed by a rectangle. The shaded region repre- 
sents a unit cell of the strip. 

To take into account the inhomogeneity introduced by 
the lost translational symmetry in the y-direction, we 
flrst express the "H^ as a sum of bond terms, i.e., as 



where //^^ (e"'^') is the spinon creation (rotor annihi- = X^h,^ ^ij- Now we treat the bonds (such as those 

lation) operator. The resulting Hamiltonian {'Hsr) in shown in Fig. 1) one at a time as our cluster and the 
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rest as the bath and obtain a mean-field two-site cluster 
Hamiltonian for each bond. Each of the bond problem 
then solved to obtain a self-consistent set of solutions 
for {^y}, y — 1, . . . , L. We treat all the unique bonds 
(coloured boxes in Fig. 1) in the super-cell (indicated by 
the shaded region in Fig. 1) utilizing the translational 
symmetry in the x-direction. The above inhomogeneous 
calculation becomes numerically expensive (for typical 
width of L = 100 considered here) as it involves diag- 
onalization of a large number of rotor Hamiltonians in 
each iteration of the self consistency loop. 

We also obtain the bulk phase diagram (lattice with 
no boundaries) where the parameters $, B, x are all site 
independent (homogeneous). For this, it suffices to solve 
the rotor mean-field Hamiltonian only for a single cluster. 
In the results discussed below all energies are measured in 
units of appropriate t both for the KM and BHZ models. 




FIG. 2. Bulk phase diagram - (a) KM model: There are three 
phases:- Correlated Topological Insulator (CTI), DSL (Dirac 
Spin Liquid), STI (Spinon Topological Insulator) phases. 
DSL/STI phases are Mott insulating, (b) BHZ model: There 
are three phases:- CTI, Correlated Band Insulator (CBI) and 
MI (Mott insulator). 

KM Model Results: We first discuss the bulk phase 
diagram of the KM model (see Fig. 2(a)) obtained within 
our formulation. For any given A there is a critical Ub 
such for U > Ub a Mott insulating state is obtained. For 
U <Ub (region below the red line Fig. 2(a)), there are lo- 
cal charge fluctuations even though there is a charge / spin 
gap, and the spinon dispersion is topologically nontriv- 
ial - we call this state a correlated topological insulator 
(CTI). There are two regimes of A for which the Hubbard 
interaction produces different types of Mott states. For 
A < Ac = |, the Mott state is a Dirac spin hquid (DSL), 
i. e., a state with no charge fluctuations and a gapless 
spinon dispersion similar to that of graphene. In this 
regime effective A for the spinon hopping renormalizes to 
zero. For A > Ac, the Mott state has a gapped topo- 
logically non-trivial spinon spectrum (spinon topological 
insulator (STI)) - the state has local spin fluctuations 
but no charge fluctuations. All the boundaries in Fig. 2 
correspond to second order transitions in that $ goes to 
zero continuously with increasing U. 

We now turn to zigzag edge terminated ribbons of finite 
width. Fig. 3 shows plot of $j, as a function of U and y for 
two values of A, e.g., 0.1 and 2. For A = 0.1 < Ac, the bulk 




FIG. 3. Mott transition in KM ribbons of width L = 50. in 
the two different regimes of A. (a) Synchronous Mott tran- 
sition obtained for A = 0.1 where the entire ribbon becomes 
Mott insulating at a critical value of Uc- The bulk value of 
critical Hubbard interaction is Ub = 6.19. (b) Asynchronous 
Mott transition for A = 2 {Ub = 17.8). The edge sites go 
Mott insulating at Ue = 15.8, followed by successive layers 
becoming Mott insulating with increasing U until the ribbon 
becomes fully Mott insulating a.t U = Uc- Dependence of Uc 
on the width L of the ribbon is shown in the upper insets in 
both cases. 

value of the critical Hubbard interaction is Ub — 6.19. As 
shown in Fig. 3(a), for this value of A, the Mott transition 
of the ribbons occurs in a synchronous'^ fashion, i. e., 
the charge fiuctuations vanish throughout the strip at 
a critical value oi Uc- Uc is quite close to Ub and the 
difference between Ub and Uc falls with increasing width 
of the ribbon (see inset in Fig. 3(a)). 

For the value of A = 2 > Ac we find completely differ- 
ent physics. In this case, the edge sites (see Fig. 3(b)) 
undergo a local Mott transition at a value of Ue — 15.8 
which is significantly smaller than the bulk value for the 
Mott transition Ub = 17.8 while the sites in the bulk 
continue to enjoy charge fluctuations. Most interestingly, 
this critical value U e does not depend on the with L of 
the ribbon over the range of ribbon widths L = 10 to 100 
studied in this work. Further increase of U above Ue 
results in successive layers attaining the Mott insulating 
state (Fig. 3(b)) - a phenomenon we call asynchronous^^ 
Mott transition. The process of asynchronous Mott tran- 
sition continues until a width dependent critical Uc is at- 
tained. Uc is always less than Ub (see inset of Fig. 3(b)) 
and Uc ^ Ub &s L becomes larger. 

The key point to be noted is that the nature of the 
Mott transition in KM ribbons is determined by A. For 
all A < Ac we find synchronous Mott transition, while 
the second regime A > Ac, asynchronous behaviour is 
obtained. 

BHZ Model Results: The results discussed here are 
for i = u = 1. Fig. 2(b) shows the bulk phase diagram 
of the BHZ model. There are again two regimes of the 
topological parameter. For eo < 3/2, there are two tran- 
sitions with increasing U - In this regime first transition 
is at Ua, where the spinon dispersion becomes topologi- 
cally trivial, but the system is still not Mott insulating, 
i. e., the system effectively is a correlated band insula- 



4 




0.5 



-0.25 





1 ^0 = 2.0 


(b) 






•0.10 - 




■0.05 - 


I 1 


-0.03 





0.01 



0.25 



0.5 



FIG. 4. (a) Asynchronous Mott transition in a BHZ ribbon 
with L = 100. (b) Wave function of mode that appears at 
the MI - CTI interface in the asynchronous Mott transition. 
For U < Ue, the mode is localized at the edge of the sample 
(red curve). With increasing U > Ue the outer layers become 
Mott insulating and the interface modes appear in the inside 
layers of the ribbon. 



tor (CBI). At a larger critical value Ub, the CBI state 
undergoes a Mott transition and obtains paramagnetic 
Mott insulator. The spinon dispersion renormalizes to 
zero. When eq > 3/2, the CTI phase gives way to the 
MI phase without the intervening CBI phase. The Mott 
transition in this case is a first order transition (unlike in 
the KM case). 

Fig. 4(a) shows asynchronous Mott transition in a BHZ 
ribbon for cq = 2. In fact, all BHZ ribbons that we have 
studied over a range of cq and widths L undergo asyn- 
chronous transitions. All other qualitative features of the 
transition are similar to that of the KM model. An in- 
teresting aspect is that in both models, the asynchronous 
Mott transition leads to a state with trivial topology. 
Thus between Ue and Uc the ribbon consists of two re- 
gions, one adjoining the edges that is Mott insulating and 
topologically trivial, while the center portion continues to 
be a CTI endowed with a non-trivial topology. Thus one 
expects interface modes to appear at the layer that sep- 
arates the two regions which is now in the bulk of the 
ribbon and to evolve further into the bulk as U increases 
from Ue to Uc- Indeed, we do find such spinon modes as 
shown for the BHZ model in Fig. 4(b). KM model also 
has similar physics in the asynchronous regime. 




FIG. 5. Local site compressibility nij/)/ kb for the KM model. 

Why are there two types of transitions and what gov- 
erns this? The nature of the transition is governed by the 
local nature of the topological edge states. The quantity 
K{y) = lim,3^oo dT{ny{T)ny{0)) - (ny(0))^, where /3 is 
inverse temperature, calculated using the non-interacting 



Hamiltonian is a measure of the local compressibility of 
a site in a layer at a given y. If the edge states brought 
about by the topological dispersion results in a large com- 
pressibility compared to the bulk, then the edge sites 
resist Mott localization - topological resilience to Mott 
localization - since much kinetic energy is lost in local- 
izing carriers at the edges. Indeed for the KM model, 
the velocity of edge states for A ^ Ac is 6A while for 
A = 2 the edge state velocity is 0.75A. This is reflected 
in a large site compressibility at the edges compared to 
the bulk value (kb) for A = 0.1 (see Fig. 5). For highly 
compressible edges, Mott transition of the ribbons occurs 
at value of U close to that of the bulk value and occurs 
in a synchronous fashion. In the other limit where the 
edge site compressibility is comparable or lower to that 
of the bulk (as is the case for A = 2 in KM, and for the 
BHZ model), Hubbard energy dominates since there is 
no significant kinetic energy to be gained and Mott lo- 
calization occurs at value of U significantly lower than the 
bulk. The kinetic energy of the bulk spinon modes comes 
to play; these modes prevent the synchronous Mott tran- 
sition and rendering it asynchronous. Our arguments, if 
applied, to the usual one band square lattice half-filled 
Hubbard model (no edge states) will predict an asyn- 
chronous transition for finite ribbons - indeed we do find 
this in our SRMFT simulations of such systems. 

Interestingly, for the synchronous Mott transition the 
spinon edge mode velocity is continuously renormalized 
to lower values as U approaches Uc ■ In the asynchronous 
case, the edge modes "switch layers" and eventually van- 
ish at Uc- This provides a clear physical picture of the 
nature of the Mott transition in finite ribbons of TIs. An- 
other important outcome of this study pertains to the use 
of surface probes to investigate the electronic state cor- 
related materials (particularly topological insulators). A 
key point to be borne in mind is that surface may show 
insulating character while the bulk may still be locally 
compressible. 
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